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ABSTRACT 

DARPA’s  OrbitOutlook  aims  to  augment  the  performance  of  the  Space  Surveillance  Network’s  Space  Situational 
Awareness  (SSA)  by  adding  more  data  more  often  from  more  diverse  sources  to  increase  Space  Domain  Awareness 
(SDA)  and  determine  when  satellites  are  at  risk.  OrbitOutlook  also  seeks  to  demonstrate  the  ability  to  rapidly 
include  new  instruments  to  alert  for  indications  and  warnings  of  space  events.  [1]  In  order  to  accomplish  this,  an 
underlying  framework  to  process  sensor  data  and  provide  useful  products,  such  as  angles-only  measurements,  is 
required.  While  existing  optical  signal  processing  implementations  are  capable  of  converting  raw  sensor  data  to 
angles-only  measurements,  they  may  be  difficult  to  customize,  obtain,  and  deploy  on  low-cost  demonstration 
programs.  The  Ground  Optical  Signal  Processing  Architecture  (GOSPA)  (inexpensive,  easy  to  customize,  obtain, 
and  deploy)  can  ingest  raw  imagery  and  telemetry  from  a  space-based  optical  sensor  and  perform  a  background 
subtraction  process  to  remove  anomalous  pixels,  interpolate  over  bad  pixels,  and  dominant  temporal  noise.  Streak 
end  points  and  target  centroids  are  located  using  the  AFRL  Generalized  Electro-Optical  Detection,  Tracking,  ID, 
and  Characterization  Application  (GEODETICA);  an  image  data  reduction  pipeline  for  ground  and  space-based 
SSA.  These  identified  streak  locations  are  fused  with  the  corresponding  telemetry  to  determine  the  Right  Ascension 
and  Declination  (RA/DEC)  measurements.  The  measurements  are  processed  through  the  AFRL  Constrained 
Admissible  Region  Multiple  Hypothesis  Filter  (CAR-MHF)  to  determine  its  initial  orbit  (IOD).  GOSPA 
performance  is  demonstrated  using  non-rate  tracking  collections  against  a  satellite  in  GEO,  simulated  from  a  visible 
optical  imaging  sensor  in  a  polar  LEO.  Stars,  noise  and  bad  pixels  are  simulated  based  on  look  angles  and  sensor 
parameters.  Under  this  scenario,  GOSPA  generated  an  IOD  with  absolute  position/velocity  errors  of  less  than  2.5  km 
and  0.3  m/s  using  sensor  RA/DEC  RSS  pointing  uncertainties  as  large  as  70.7  arcsec.  This  demonstrates  the 
capability  for  GOSPA  to  support  contributing  SSA  space-based  optical  sensors  under  OrbitOutlook  and 
subsequently  provide  useful  data  back  to  space  command  and  control. 

1.  INTRODUCTION 

OrbitOutlook  (O2)  is  a  DARPA  sponsored  program  with  the  objective  of  augmenting  the  Space  Surveillance 
Network  (SSN)  by  adding  contributing  optical  and  radar  sensors  with  diverse  geographic  locations.  [1]  This  includes 
space-based  as  well  as  ground-based  measurement  collection  capabilities.  While  SpaceView  exists  to  provide  a 
mechanism  to  allow  for  O2  to  remotely  control  ground  telescopes,  a  cost  effective  rapidly  deployable  signal 
processing  architecture  has  not  been  identified  or  socialized  within  the  SSA  community  to  allow  for  the  performance 
prediction  modeling  and  real-time  operational  processing  of  Resident  Space  Object  (RSO)/target  collections  using 
contributing  space-based  optical  sensors. 

While  contributing  SSA  ground-based  optical  sensor  collectors  have  a  low  cost  of  acquisition,  most  systems  are 
sensitive  to  the  background  sky  brightness  and  are  relegated  to  night  time  operation.  Space-based  optical  sensors, 
with  lower  Sun  and  Earth  exclusion  limits  can  help  to  eliminate  this  daytime  collection  gap  imposed  by  ground- 
based  optical  systems  and  are  therefore  critical  to  achieving  better  response  and  warning  times  to  space  events. 
Studies  from  the  Space-Based  Visible  program  indicate  that  with  the  use  of  efficient  search  patterns,  it  is  possible 
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demonstrate  the  ability  to  rapidly  include  new  instruments  to  alert  for  indications  and  warnings  of  space 
events.  [1]  In  order  to  accomplish  this,  an  underlying  framework  to  process  sensor  data  and  provide  useful 
products,  such  as  angles-only  measurements,  is  required.  While  existing  optical  signal  processing 
implementations  are  capable  of  converting  raw  sensor  data  to  angles-only  measurements,  they  may  be 
difficult  to  customize,  obtain,  and  deploy  on  low-cost  demonstration  programs.  The  Ground  Optical  Signal 
Processing  Architecture  (GOSPA)  (inexpensive,  easy  to  customize,  obtain,  and  deploy)  can  ingest  raw 
imagery  and  telemetry  from  a  space-based  optical  sensor  and  perform  a  background  subtraction  process  to 
remove  anomalous  pixels,  interpolate  over  bad  pixels,  and  dominant  temporal  noise.  Streak  end  points  and 
target  centroids  are  located  using  the  AFRL  Generalized  Electro-Optical  Detection,  Tracking,  ID,  and 
Characterization  Application  (GEODETIC A);  an  image  data  reduction  pipeline  for  ground  and 
space-based  SSA.  These  identified  streak  locations  are  fused  with  the  corresponding  telemetry  to  determine 
the  Right  Ascension  and  Declination  (RA/DEC)  measurements.  The  measurements  are  processed  through 
the  AFRL  Constrained  Admissible  Region  Multiple  Hypothesis  Filter  (CAR-MHF)  to  determine  its  initial 
orbit  (IOD).  GOSPA  performance  is  demonstrated  using  non-rate  tracking  collections  against  a  satellite  in 
GEO,  simulated  from  a  visible  optical  imaging  sensor  in  a  polar  LEO.  Stars,  noise  and  bad  pixels  are 
simulated  based  on  look  angles  and  sensor  parameters.  Under  this  scenario,  GOSPA  generated  an  IOD 
with  absolute  position/velocity  errors  of  less  than  2.5  km  and  0.3  m/s  using  sensor  RA/DEC  RSS  pointing 
uncertainties  as  large  as  70.7  arcsec.  This  demonstrates  the  capability  for  GOSPA  to  support  contributing 
SSA  space-based  optical  sensors  under  OrbitOutlook  and  subsequently  provide  useful  data  back  to  space 
command  and  control. 
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for  contributing  SSA  sensors  (especially  slew  rate  limited  payloads)  to  cover  most  of  the  GEO  belt  within  six  hours. 
[3]  Subsequently,  with  the  enhanced  SSA  collection  capabilities  offered  by  O2,  the  space  command  and  control 
would  not  be  able  to  process  the  increased  amount  of  data  that  could  be  delivered  and  it  could  not  be  readily  used.  [2] 
In  order  to  obtain  useful  information  from  an  optical  sensor  collection,  the  target  must  be  distinguished  from  noise 
and  other  signal  sources  (e.g.  stars,  planets,  and  other  space  objects).  This  type  of  processing  involves  detailed 
knowledge  of  the  optical  detector  as  well  as  the  tracking  CONOPS,  platform  position  and  orientation  during  the 
collection  period,  and  the  ability  to  associate  objects  obtained  from  separate  observations. 

The  Millennium  Space  Systems  led  Ground  Optical  Signal  Processing  Architecture  (GOSPA)  addresses  the  current 
need  for  a  low-cost  Space  Situational  Awareness  (SSA)  processing  architecture  solution  in  support  of  DARPA 
TTO’s  O2  program  to  include  systems  that  were  not  originally  designed  to  contribute  SSA  but  are  technically 
capable.  Fig  1  summarizes  the  physics  based  simulation  which  was  developed  for  the  verification  and  validation  of 
the  GOSPA  process. 


Fig  1:  Physics  Based  Optical  Models  and  High-Precision  Orbit  Propagation  Integrated  into  the  GOSPA  Simulation 


2.  SPACE-BASED  OPTICAL  SENSOR  MODEL 


A  hypothetical  contributing  SSA  sensor  using  a  generic  representation  of  a  body  fixed  space-based  sensor  is 
proposed  [3].  The  optical  parameters  and  sensor  properties  used  in  the  Signal  to  Noise  Ratio  (SNR)  model  are 
provided  in  Table  1. 

Table  1:  Optical  Sensor  Modeling  Parameters  with  accompanying  graph  of  sensor  performance  for  a  space-based 
LEO  sensor  observing  a  GEO  object  assuming  a  relative  angular  rate  of  30  arcsec. 

Minimum  Detectable  RSO  vs  Solar  Phase  Angle 

290 


240 


190 


140 


90 


40 

0  50  100 


Focal  Plane  Array  Size 

512  x  512  pixels 

Pixel  Size  (Length  x  Width) 

20  pm  x  20  pm 

Well  Capacity  (electrons) 

170,000 

Digital  Quantization  (A/D  bits) 

12  bits 

Quantum  Efficiency 

80% 

Optical  Throughput 

0.80 

Integration  Time 

1  sec 

Spectral  Range 

0.300  -  0.900  pm 

Sensor  Field  of  View 

1.0  x  1.0  deg 

Frame  Stacking 

1  Co-added  Frames 

Aperture  Diameter 

15  cm 

Aperture  Area  Obscuration 

0.10 

Number  of  Frames  per  Frameset 

8  Frames 

Instantaneous  Field  of  View 

7.04  arcsec/pixel 

Energy  on  Detector  (EOD) 

0.70 

Read  Nosie 

50  electrons 

Pixel  Area  Fill  Factor 

1 

Gain 

12 

Reset  or  Calibration  Time 

20  seconds 

Multiplexer  Voltage  (full  swing) 

2  volts 

Dark  Current 

7E- 1 8  Amps 

Dark  Current/Responsivity  Non- 
Uniformity 

0.00002  Fraction 

Background  Visual  Magnitude 

22  Mv/arcsec2 

Probability  of  False  Alarm  (Pfa) 

IE-4 
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The  configuration  of  the  optical  sensor  model  is  illustrated  in  Fig  2. 


Fig  2:  Optical  Sensor  Model 


The  sensor  model  has  the  capability  of  converting  visual  magnitudes  to  radiometric  quantities.  However,  since  the 
model  covers  the  full  range  of  light  (from  UV  through  VLWIR)  a  more  generic  description  of  this  model  is 
provided.  An  optical  signal,  due  to  a  point  source  or  streak,  is  given  by  (1) 


signal  =  ®tot  •  Aoptlcs  •  rj^  ■  EODeffective  •  QE  •  tml 


(1) 


Where  ®  is  the  total  photon  flux  of  the  target  at  the  aperture  of  the  optical  system,  Aoptics  is  the  aperture  area. 


Coptics ’s  the  average  optical  transmission  per  optical  element,  N.  EODeffective  is  the  effective  energy  on  detector  for 


a  streaked  or  stationary  target,  QE  is  the  in-band  quantum  efficiency,  and  tint  is  the  integration  time  of  each  frame. 
The  signal  to  noise  ratio,  SNR,  is  given  by  (2) 


SNR  = 


signal 

noise 


coadds 


(2) 


Where  Ncoadds  is  either  the  number  of  frame  co-adds  or  time  delay  integrations.  The  noise  is  given  by  (3)[4][5]: 
noise  =  ^dNphotons+ dNphotonsFPN+  dN 2i±  H-dN^™  +  dNj  +dN^  +  dN3lf  +dN;tc  +dN2vlf  +  dN^  +dN^  (3) 


Where  dNphotons  =  ^signal  +dNfheiiml  +dNfodmcal  and  dNread  is  the  read  noise  in  noise-electrons.  dNthermid  is 

the  photoelectron  noise  due  to  the  temperature  of  the  focal  plane.  dNzodiacal  is  the  photoelectrons  due  to  background 
light.  dNphotonFPN  is  the  sensor  responsivity  fixed  pattern  noise.  dNdrk  is  the  detector  dark  current  noise.  dNdrkFPN  is  the 
detector  dark  current  non-uniformity.  dNj  is  the  Johnson  and  Shunt  Noise.  dNdlf  is  the  detector  1/f  noise.  [6]  dNktc  is 
the  capacitance  switching  noise.  dNvlf  is  the  amplifier  input  referred  1/f  noise.  [5]  dNaw  is  the  amplifier  input 
referred  white  noise,  and  dNa2d  is  the  digitization  noise.  Note  that  unless  otherwise  specified,  the  units  of  all 
individual  noise  terms  are  in  electrons  or  noise-electrons.  Other  functions  include  an  approximation  for  the 
resistance-area  product,  R0A,  vs  detector  cutoff  wavelength,  the  effective  EOD,  and  the  formulas  for  the  various 
radiometric/photometric  quantities. 


3.  ORBIT  PROPAGATION  FORCE  MODEL 


The  General  Mission  Analysis  Tool  (GMAT),  developed  by  the  Goddard  Space  Flight  Center  at  NASA,  was  chosen 
to  perform  the  orbit  propagation  of  both  the  target  satellites  and  the  sensor  platform.  For  this  analysis,  it  is  assumed 
that  the  targets  and  the  sensor  are  not  maneuvering  during  the  simulation  period.  The  second  order  differential 
equation  to  describe  the  motion  of  a  satellite  is  provided  (4)  [7]: 


d2r 

d2t 


X 


m 


~4  r+Vfcj 

r 


+ 


cl 


m, 


k=  1 
k*j 


JL 

r3 


p  r  A 

_l_  _  sR^R-^e  r 
772 


sG 


(4) 
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Where  —  r  describes  the  central  body  point  mass,  V  (j).  is  the  central  body  direct  non-spherical  force  component. 


*= i 

k*J 


,3  ,-3 

V  r*/  J 


;  the  direct  and  indirect  third  body  point  mass,  and  PsrCr\>  p  describes  the  solar  radiation 


pressure. 

For  simulation  purposes,  the  following  force  model  parameters  and  models  were  chosen  for  the  propagation  of  both 
sensor  and  target  satellites: 


Table  2:  General  Mission  Analysis  Tool  Force  Model  Settings  for  All  Satellites 


Central  Body  Point  Mass 

Earth 

Third  Body  Point  Mass 

Sun,  Mercury,  Venus,  Moon,  Mars,  Jupiter, 

Saturn,  Uranus,  Neptune,  and  Pluto 

Earth  Gravity  Model 

Earth  Gravity  Model  1996 

Gravity  Field  Degree  and  Order 

100  x 100 

Spacecraft  Coefficient  of  Reflectivity  (Cr) 

0.7 

f  107 Average  (81 -day  averaged  10.7cm  solar  flux) 

150 

F107daily  (Daily  10.7cm  solar  flux) 

150 

Daily  magnetic  flux  index 

4 

Solar  Radiation  Pressure  Area  (A) 

10  m2 

Spacecraft  mass  (m) 

1000  kg 

For  this  simulation,  the  space-based  sensor  platform  was  given  a  polar  Low  Earth  Orbit  described  by  the  initial  ECI 
state  vector  provided  in  Table  3. 

Table  3:  Space-Based  Sensor  Initial  ECI  State  Vector 
Epoch  Date  (UTC)  x-position  y-position  z-position  x-velocity  y-velocity  z-velocity 

(km)  (km)  (km)  (km/s)  (km/s)  (km/s) 

I  01/01/2014  00:00:00  I  7578.13700  I  0  T  0  I  0  I  0  I  7.2525  ] 


The  initial  ECI  state  vectors  of  the  Geosynchronous  Earth  Orbit  (GEO)  targets  in  which  the  sensor  is  tasked  to 
obtain  tracks  are  included  in  Table  4. 


Table  4:  GEO  Target  ECI  State  Vectors 


Epoch  Date  (UTC) 

x-position 

(km) 

y-position 

(km) 

z-position 

(km) 

x-velocity 

(km/s) 

y-velocity 

(km/s) 

z-velocity 

(km/s) 

01/01/2014  00:00:00 

42164.00000 

0 

0 

0 

3.07467 

0 

01/01/2014  00:00:00 

42163.98395 

36.79503 

0 

-0.00268 

3.07467 

0 

01/01/2014  00:00:00 

42163.98395 

-36.79503 

-36.79501 

0.00268 

3.07467 

0 

4.  TARGET  ILLUMINATION  MODEL 


For  model  simplification,  all  targets  are  modeled  as  a  Lambertian  sphere.  The  total  photon  flux  (i.e.  illumination)  of 
the  target  is  dependent  on  the  waveband  of  the  observing  sensor  in  which  to  calculate.  For  the  generic  space-based 
optical  sensor,  the  waveband  is  in  the  visible  spectrum,  note  that  this  target  illumination  model  will  support  other 
wavebands  as  well.  Also,  because  the  target  is  in  a  Geosynchronous  Orbit,  its  distance  from  Earth’s  is  large  enough 
such  that  the  Earth  can  be  modeled  as  a  Lambertian  sphere  as  well.  This  simplifies  the  Earth  reflectance  models  and 
allows  the  use  of  a  Lambertian  spherical  approximation. 

The  total  photon  flux,  ®tot,  in  photons  s"1  m'2  of  a  target  observed  from  a  sensor  can  be  found  in  (5): 


tot 


^TSE  ®ES  ^TE 


(5) 


Where  ®tse  is  the  photon  flux  from  self-emission  (thermal)  of  the  target,  ®es  is  the  photon  flux  from  the  Earth 
shine  reflecting  on  the  target,  ®s  is  the  photon  flux  from  the  solar  reflection  on  the  target,  and  ®te  is  the  photon 
flux  from  the  thermal  radiation  reflected  from  the  Earth  on  the  target. 

The  photon  emissivity  (P  s'1  m'2)  from  thermal  self-emission  given  a  specific  sensor  waveband  can  be  numerically 
computed  from  the  Planck  function  (6)  [8], 
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°tse  L 


exp(c2/lJ '  )-l 


Where  c ,  and  c2  are  constants,  c,  =  2 zdic2  =  3.74 x  10  16  Wm  2  and  c  =  _  =  \  44xl0~2m  K 

2  k 

The  temperature,  Ttgt,  of  the  target  can  be  approximated  using  an  isothermal  uniform  satellite  temperature  model 
using  the  conservation  of  energy  equation  (7)  (i.e.  energy  emitted  is  equal  to  the  energy  absorbed)  [9]: 

£tgt°AtgtTtgt  =  £tg1°AtgtFtgt,eTe  +  Qsun  +  Qer  +  Qi  ^ 

Rearranging  (7)  for  the  uniform  temperature  of  the  target,  in  steady  state,  is  found  from  (8): 

I  a  Q  +Q  +  0 

Ttgt  =  4  Ftgt,eTe4  + Vsun  Ver  Vl  (8) 

V  stgtaAtgt 

Where  £lgtis  the  target  emissivity.  G  is  the  Stefan-Boltzman  constant.  Teis  the  uniform  temperature  of  the  Earth.  Q, 

is  the  internally  generated  power  of  the  target  in  watts  (from  averaging  the  generated  power  of  many  active  satellites 
over  their  spherical  area  in  GEO,  assume  271  W/m2).  Qer  is  the  Earth  reflected  solar  input  which  can  be  computed 
analytically  for  satellites  in  GEO  as  shown  in  equation  (9): 

Qer  =  I»unaE  (“tgt  A±  )  (9) 

TTClj 

The  phase  integral  function,  P(0),  can  be  described  in  (10): 

p(9)  =  ((jt  -  0)cos0  +  Sin9)  (io) 

371 

A_  is  the  target  cross  section  area,  OCtgt  is  its  surface  absorptivity,  Isun  ,  is  the  total  solar  intensity  value  of  the  sun 
reaching  Earth’s  surface  (1366.1  W/m2  corresponds  to  the  STME  E490  2000  Standard),  di  is  the  distance  from  the 
center  of  the  earth  to  the  target,  AE  is  the  surface  area  of  the  Earth,  (XE  is  the  Earth  Albedo,  and  0X  is  the  angle 
between  the  Earth-Sun  and  the  Earth-target  vectors.  The  solar  input  to  the  spacecraft  is  described  by  (11). 

Qsun  =  ttlgtAdsun  (U) 

Ftgt.e  is  the  view  factor  of  the  Earth  from  the  target  and  can  be  computed  in  (12). 


R2  1 

=  1-  1 - f  - 


The  photon  flux  from  reflected  solar  energy  from  the  Earth  is  described  by  (13) 


^ES  ^Asun 


aEAEp(0,)atgtA±p(02) 


Where  1  Asun  is  the  irradiance  of  the  sun  given  a  specific  waveband  (can  be  determined  from  a  table  lookup  given 

from  the  2000  ASTM  Standard  Extraterrestrial  Spectrum).  d2  is  the  distance  from  the  target  to  the  sensor.  02  is  the 
angle  between  the  two  vectors  that  describe  d  and  d2  (Earth  phase  angle  on  the  target).  The  flux  from  direct  solar 
illumination  can  be  determined  from  (14): 

attALp(<f) 

=  >.sun - —2 -  (14) 


Where  E,  is  the  solar  phase  angle  (angle  between  the  sun-target  vector  and  the  target-sensor  vector,  d2).  The  photon 
flux  related  to  the  thermal  spectrum  reflectance  of  the  Earth  temperature  on  the  target  can  be  found  from  (15): 
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Earth  Thermal 


(15) 


(Dte  —  I 


«,gtA±P(92) 
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Where  IgaithThennai 's  the  computed  flux  over  the  sensor  waveband  using  the  thermal  self-emission  in  (6). 

5.  SENSOR-TARGET  TRACKING  CONOPS 


The  space-based  platform  must  slew  to  align  the  sensor  boresight  vector  in  the  direction  of  the  object  it  has  been 
tasked  to  track.  Since  it  is  assumed  that  there  is  no  additional  tracking  intelligence  on-board  the  sensor  payload, 
consider  the  scenario  where  the  sensor  is  commanded  to  slew  to  a  certain  RA/DEC  position  at  a  specified  time.  This 
is  considered  a  no-track  CONOPS  as  opposed  to  a  sidereal  track  CONOPS  (e.g.  where  the  sensor  is  tracking  the 
relative  star  motion)  or  an  object  track  CONOPS  (e.g.  where  the  sensor  is  tracking  the  relative  motion  of  the  object). 
This  no-track  CONOPS  appears  to  be  representative  of  the  type  of  tracking  performance  one  would  expect  when 
working  with  space-based  sensors  designed  for  other  missions  outside  of  SSA.  Tracking  opportunities  were  selected 
based  on  predicting  the  SNR  and  Probability  of  detection,  Pd,  of  a  target  of  interest  and  selecting  the  opportunities 
corresponding  to  the  greatest  Pd  without  violating  any  sensor  constraints  (i.e.  sensor  solar  exclusion  angles, 
maximum  slew  rate  constraints,  etc.) 

6.  SPATIAL/TEMPORAL  BACKGROUND  AND  ANOMALOUS  PIXEL  CORRECTION 


Under  certain  circumstances,  the  data  obtained  from  an  electro-optical  sensor  system  has  sufficient  noise  such  that 
pre-processing  is  necessary  before  running  a  feature  point  detection  algorithm  on  the  collection  data.  Under  those 
circumstances,  further  processing  is  performed  to  reduce  random  noise  (e.g.  frame  stacking,  velocity  filtering). 

Some  performance  tuning  may  be  necessary  as  this  has  an  associated  artifact  of  further  reducing  the  signal  strength. 
At  a  very  high  level,  it  is  possible  to  combine  collection  data  from  different  points  in  time  and  slightly  different  pixel 
positions  together  and  use  that  knowledge  to  further  clean  each  individual  frame.  This  process  is  outlined  in  Fig  3. 


Temporally/Spatially  C 

ffset  Reference  Frames 

Temporally/Spatially  1 

Iffset  ReferenceFrames 

1 

/ 

/ 

/ 

/ 

/ 

/ 

- 1 

Unprocessed  Observation  Frames 


Compute  average  of  the  two  reference  frame  stacks  which 
temporally  bracket  target  frame  stack  and  subtract  from  each 
target  frame 

a.  This  is  the  0th  order  background  subtraction 

i.  Target  and  reference  frames  are  obtained 
with  the  exact  same  integration  settings  so 
dark  current,  read-noise  etc  ...  are  the  same 
Identify  spatially  anomalous  pixels  using  the  average  of  two 
reference  frame  stacks 

a.  Apply  complementary  3x3  median  spatial  filter  to 
difference  of  the  two  average  reference  frames 
Compute  robust  statistics  and  identify  exceedances 
Eliminate  exceedances  and  recomputed  robust  statistics 
to  identify  further  exceedances 
Generate  bad  pixel  map 


b. 


2)  Identify  temporally  anomalous  pixels  using  each  reference 

frame  stack  separately 

a.  Treat  each  reference  pixel  as  a  temporal  time  series 

b.  Compute  temporal  covariance  matrix  averaged  over  all 
pixel  time  series 

c.  Use  eigenvector  decomposition  of  covariance  matrix  to 
compute  Mahalanobis  distance  for  each  time  series 

i.  Mahalanobis  distance  is  the  N-dimensional 
expression  of  an  exceedance  in  units  of  o 
accounting  for  covariance 

d.  Reject  pixels  with  Mahalanobis  distance  exceeding 
specified  threshold  (i.e.  4o) 

e.  Repeat  from  step  (b)  excluding  outlier  pixels  until 
convergence  criteria  met  (e.g.  covariance  matrix  is 
stable  within  a  specified  tolerance) 

f.  Generate  a  bad  pixel  map 


Use  Delaunay  triangulation  to  interpolate  over  bad  pixels 
identified  above  for  each  background-subtracted  target  frame 
Simple  spatial  filter  using  target  frame  stack 

a.  Assuming  optical  blur  smears  target  signal  over  several 
pixels,  apply  complementary’  3x3  median  spatial  filter 

b.  Compute  frame  statistics  and  identify  exceedances  as 
bad  pixels 

c.  Replace  bad  pixels  with  median  of  3  x  3  neighborhood 


Eliminate  correlated  temporal  noise  in  target  frame  stack: 

a.  Treat  each  target  pixel  as  a  temporal  time  series  4) 

b.  Compute  temporal  covariance  matrix  averaged  over  all 

pixel  time  series  6) 

c.  Compute  eigenvector  decomposition  of  covariance 
matrix 

i.  Eigenvector  with  the  largest  eigenvalue 

contains  majority  of  power  associated  with 
common  mode  temporal  modulation 

d.  Eliminate  dominant  eigenvector  and  re-project  pixel 
time  series  into  reduced  subspace  to  eliminate  common 
mode  signal 

i.  Assuming  equipartition  of  target  signal 

among  eigenmodes,  this  will  reduce  target 
signal  by  1/total  number  of  eigenmodes,  but 
in  practice  will  improve  SNR  by  reducing 
the  noise  by  an  even  larger  factor.  Pre-Processed  Observation  Frames 

Fig  3:  Spatial/Temporal  Background  &  Anomalous  Pixel  Correction  Process 
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7.  GEODETICA:  RSO  DETECTION  AND  DISCRIMINATION  ALGORITHM 


Once  a  set  of  frames  has  been  suitably  pre-processed,  space  objects  are  detected  and  localized  within  each  frame 
using  a  segmentation  algorithm  based  on  the  concept  of  Phase  Congruency  followed  by  a  centroid  and  corner 
localization  procedure.  In  image  processing  segmentation  is  the  assignment  of  image  pixels  as  belonging  to  either 
foreground  objects  or  the  background.  Often  in  the  processing  of  images  from  space  sensors  image  segmentations 
are  performed  by  setting  a  threshold  in  image  intensity.  Global  thresholds  in  intensity  cause  many  undesirable 
effects  due  to  non-flat  noise  backgrounds  caused  by  stray  light  or  other  artifacts  from  pre-processing  which 
complicate  subsequent  RSO  localization  and  discrimination  processes.  Other  algorithms  which  perform  this 
function  also  perform  segmentation  on  various  additional  aspects  of  image  data  including  segmentation  based  on 
motion,  signal  to  noise  ratio,  and  higher  order  statistical  moments.  The  segmentation  based  on  Phase  Congruency 
used  in  this  paper  assigns  pixels  based  on  moments  of  a  measure  of  local  energy  which  is  invariant  to  local  image 
contrast.  Additionally,  subsequent  localization  of  segmented  signals  do  not  suffer  the  biases  typically  caused  by 
other  local  energy  segmentation  methods  due  to  the  implicit  Gaussian  smoothing  used.  The  Phase  Congruency 
model  for  feature  detection  does  not  assume  features  are  local  maxima  in  the  intensity  gradient,  rather  it  postulates 
that  features,  in  this  case  star  and  RSO  dots  and  streaks,  are  located  at  points  in  the  image  where  the  Fourier 
components  are  maximally  in  phase.  The  local  energy  of  a  signal,  PC(x),  is  represented  using  (16): 


pcOO  = 


|E(x)| 


(16) 


£nAn(x) 

The  standard  way  of  visually  conveying  this  concept  is  to  notice  that  the  Fourier  components  are  all  in  phase  at  the 
location  of  the  step  edge  in  the  function  shown  in  Fig  4  shows  an  additional  geometric  interpretation  of  the  concept 
where  if  the  Fourier  components  are  plotted  tip  to  tail  as  complex  vectors,  an  in  phase  and  therefore  perceptible 
location  would  achieve  a  large  radius  or  have  a  large  local  energy,  while  a  location  which  is  more  like  noise  would 
behave  like  a  random  walk  and  not  achieve  the  same  value  of  local  energy.  In  fact  this  concept  can  be  used  to 
model  the  behavior  of  noise  in  a  given  sensor  and  to  define  thresholds  which  are  directly  related  to  the  probability  of 
false  alarms. 

Imaginary 


Fig  4  (left)  Multiple  Fourier  Components  of  a  periodic  signal  in  phase  at  the  step  edge,  (right)  Fourier  components 

plotted  tip  to  tail  in  the  polar  diagram 


The  concept  of  Phase  Congruency  and  its  application  to  imagery  for  the  purposes  of  finding  edges  and  corners  was 
first  developed  by  Peter  Kovesi  [10].  In  that  work  he  explains  that  the  minimum  and  maximum  moments  of  this 
measure  can  be  directly  computed  and  used  to  define  edges  and  corners.  For  the  work  performed  in  this  paper,  the 
corners  are  used  to  localize  space  object  centers  which  are  typically  co-located  with  intensity  centroids  of  constant 
intensity  signals,  as  well  as  streak  endpoints  in  the  case  that  apparent  motion  was  captured  over  a  single  image 
exposure  time.  Currently  only  centroids  are  associated  frame  to  frame  but  sufficient  information  is  provided  by  the 
corners  for  additional  feature  descriptions  which  can  be  used  for  object  and  motion  discrimination.  Figure  6  shows 
the  application  of  this  segmentation  technique  followed  by  centroiding  and  corner  localization  in  order  to  extract  a 
streaking  signal  from  an  image  with  a  gradient  background.  Following  the  extraction  of  star  and  RSO  signals  in 
each  frame,  a  nearest  neighbor  association  is  performed  in  order  to  build  the  observed  track  of  the  RSO  and 
determine  the  observed  right  ascension  and  declination  as  a  function  of  time.  Once  these  parameters  have  been 
determined  for  all  observations,  an  orbit  for  the  observed  RSO  can  be  estimated. 
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Fig  5:  Streaking  RSO  Segmented  Using  GEODETICA  Including  Simultaneous  Determination  of  the  Centroid  and 

Endpoints  (Note  this  graphic  is  simulated  sensor  data) 

8.  CAR-MHF  ORBIT  DETERMINATION  ALGORITHM 


The  Constrained  Admissible  Region  Multiple  Hypothesis  Filter  (CAR-MHF)  is  an  optical  measurement  (i.e.  angles- 
only)  initial  orbit  determination  (IOD)  strategy  that  combines  the  track  initialization  capability  of  the  CAR  (which 
creates  a  set  of  hypothetical  trajectory  estimates)  with  an  MHF  that  implements  an  unscented  Kalman  filter  (UKF), 
see  Stauch  et  al.  for  a  more  detailed  discussion  [11].  As  the  name  suggests,  the  CAR  involves  constraining  the 
possible  ambiguity  space  to  an  admissible  region,  following  the  method  presented  by  Milani  [12].  A  tracklet  of 
optical  data  (a  set  of  3  or  more  closely-spaced  angles-only  measurements)  can  be  compressed  to  form  an  angles  and 
angle-rates  measurement  quadruplet  (based  on  the  raw  measurement  noise).  As  shown  by  DeMars  et  al.,  [13]  this 
measurement  can  be  used  to  create  semi-major  axis  (SMA)  and  eccentricity  (e)  constraints  in  range/range -rate  space, 
creating  a  CAR  tailored  to  the  type  of  RSO  one  might  be  interested  in  (i.e.  GEO).  Even  when  constrained,  the  region 
is  likely  too  large  to  be  used  to  initiate  a  single  filter,  thus  a  grid  can  be  formed  within  the  CAR  to  form  a  set  of 
suitable  hypotheses.  This  yields  a  set  of  hypothesis  in  observation  space,  which  can  be  passed  through  an  unscented 
transform  [14]  to  produce  a  set  of  hypotheses  in  Cartesian  state  space  that  can  readily  form  the  a  priori  means  and 
covariances  for  the  MHF. 

The  MHF  used  in  CAR-MHF  follows  the  method  introduced  by  DeMars  and  lah  [15]  in  which  the  estimate  of  an 
object  is  represented  by  multiple  hypotheses.  It  utilizes  the  UKF,  a  special  case  of  the  sigma -point  Kalman  filter 
(SPKF)  [16]  that  employs  the  unscented  Transform  [14],  As  opposed  to  the  linear  traditional  Kalman  filter,  the  full 
non-linear  equations  of  motion  and  measurement-state  equations  are  used.  Each  hypothesis  (initialized  via  CAR)  is 
propagated  and  updated  with  a  UKF,  in  parallel.  Additionally,  a  weighting  term  is  carried  with  each  hypothesis.  All 
hypotheses  are  initially  weighted  equally  and  the  weights  are  held  constant  during  the  propagation  step.  During  the 
update  step,  each  hypothesis’  weight  is  updated  based  on  the  likelihood  that  hypothesis  originated  the  observations. 
Hypotheses  may  be  pruned  if  and  when  its  weight  falls  below  some  desired  threshold,  allowing  the  MHF  to 
eventually  converge  to  a  single  hypothesis. 

9.  RSO  DETECTION  PERFORMANCE  VS  SNR  FOR  NO-TRACK  COLLECTIONS 


To  date  it  has  been  observed  that  this  architecture  can  be  used  for  the  detection  of  both  single  RSOs  as  well  as 
multiple  RSOs  within  a  collection.  Parametric  analysis  to  vary  the  SNR  of  the  RSO  was  performed  in  order  to 
bound  the  performance  of  the  GOSPA  process.  Simulation  results  indicate  that  an  RSO  with  an  SNR  above  2  can  be 
effectively  detected  and  discriminated. 


Tgt  SNR  =  11  3925  Tgt  SNR  =  3.7563  Tgt  SNR  =  2.237 
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Fig  6:  Single  RSO  Detection  Performance  (Note  this  graphic  is  simulated  sensor  data) 

Empirically,  it  is  easier  to  automatically  identify  and  track  a  single  RSO  using  the  current  implementation  of  the 
Feature  Point  Detection  Algorithm  at  low  SNR  values  than  it  is  to  identify  and  track  multiple  RSOs  at  the  same 
SNR.  It  was  found  that  multiple  RSO  detections  are  possible  assuming  an  SNR  of  6  or  above.  The  performance  of 
the  detection  and  discrimination  architecture  is  critically  associated  with  the  sensor  sensitivity,  field  of  view,  and 
other  parameters  as  the  algorithms  are  written  to  be  capable  of  discerning  if  the  sensor  is  in  a  rate  track,  sidereal 
track,  or  other  modes  which  require  a  minimum  number  of  stars  per  frame.  In  cases  where  the  number  of  RSOs  is 
larger  than  the  observed  stars,  it  can  be  difficult  to  distinguish  without  a  catalog  this  scenario  from  a  rate  track. 
Future  work  is  focusing  on  enhancing  this  capability. 

10.  CAR-MHF  PERFORMANCE  VS  POINTING  ANGLE  UNCERTAINTY 

Many  low-cost  space-based  sensor  solutions  have  much  larger  error  in  pointing  knowledge  due  to  losing  lock  on 
stars  with  a  star  tracker,  IMU  drift,  and  other  GN&C  limitations.  Without  the  use  of  star  calibration  to  enhance 
pointing  knowledge,  it  is  still  possible  to  determine  a  coarse  initial  orbit  of  a  series  of  objects  using  CAR-MHF.  For 
the  purposes  of  this  study,  CAR-MHF  was  used  to  perform  an  Initial  Orbit  Determination  from  a  single  space-based 
sensor  platform  (as  described  in  the  previous  sections).  The  RSS  positionYvelocity  3-0  uncertainty  associated  with 
each  variation  of  collection  cadence  was  determined  for  various  sensor  pointing  uncertainties.  This  should 
effectively  quantify  the  operational  performance  of  the  hypothetical  sensor  system. 

Fig  7  shows  the  root  sum  square  (RSS)  of  the  position  and  velocity  IOD  state  vector  subtracted  from  the  simulation 
state  vector  (truth).  Fig  8  displays  the  ratio  of  the  absolute  error  in  IOD  over  the  1-CT  uncertainty  estimates  of  CAR- 
MHF.  The  sensor  pointing  uncertainty  was  varied  (assuming  a  Gaussian  noise  distribution)  from  20  -  50  arcsec  in 
both  apparent  declination  and  right  ascension  (28.3-70.7  arcsec  RSS).  Results  show  that  while  it  is  possible  to  obtain 
an  IOD  on  five  collections  from  a  single  space-based  sensor  over  a  24-hour  period,  the  RA/DEC  pointing 
uncertainty  must  be  at  or  below  30  arcsec.  When  the  number  of  collections  from  a  single  space-based  sensor  of  the 
object  is  sufficiently  large  (e.g.  10-12  collections/day)  there  is  a  point  of  diminishing  returns  for  both  the  position 
and  velocity  error  of  the  IOD  as  both  states  converge  to  the  less  than  the  RA/DEC  pointing  uncertainty. 
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Fig  7:  Absolute  Position  and  Velocity  Error  (Using  Truth)  of  IOD  with  Variation  in  Pointing  Uncertainty  and 
Number  of  Collections.  3-cr  Position  Uncertainty  for  the  50  arcsec  (70.7  arcsec  RSS)  case  overlaid 
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Fig  8:  Absolute  Error  Expressed  as  #-CT  Uncertainty  of  IOD  with  Variation  in  Sensor  Pointing  Uncertainty  and  Total 

Number  of  Collections 


Note  that  for  all  cases  where  the  sensor  noise  was  above  20  arcsec,  the  absolute  error  in  the  state  vector  estimate 
falls  between  0  and  3-cr  of  the  CAR-MHF  reported  uncertainty  estimates.  And  for  all  sensor  noise  cases,  the 
absolute  position  error  is  well  below  the  reported  3-cr  uncertainty  for  the  50  arcsec  sensor  noise  pointing  case. 
Possible  causes  of  the  relatively  high  error  to  uncertainty  ratios  (see  Fig  8)  for  lower  sensor  noise  cases  are  because 
of  systematic  errors  induced  by  the  process  of  extracting  the  RSO  position  from  the  centroiding  algorithm 
(GEODETICA  is  centroiding  on  the  middle  of  each  streak  instead  of  using  one  of  its  end  points),  simulated 
discretization  of  the  target  on  a  7  arcsec/pixel  sensor  focal  plane,  different  orbit  high-precision  propagation  models 
(CAR-MHF  uses  an  internal  force  model  with  similar  but  not  identical  settings  to  GMAT),  and  the  range/range-rate 
constraints  within  CAR  do  not  yet  account  for  angle/angle -rate  noise. 

11.  CONCLUSIONS  AND  FUTURE  WORK 


This  research  has  demonstrated  that,  with  GOSPA,  it  is  possible  to  simulate  a  day  in  the  life  of  an  OrbitOutlook 
scenario  and  assess  the  IOD  performance  of  a  sensor  system  and  it’s  CONOPS  of  looking  at  a  GEO  object  from 
LEO.  An  IOD  was  obtained  with  as  few  as  four  collections  over  a  24-hour  period.  Sensor  pointing  uncertainties  as 
large  as  70.7  arcsec  RSS  were  successfully  processed  using  CAR-MHF.  Single  RSOs  with  SNRs  above  two  were 
properly  detected  using  GEODETICA.  GOSPA  is  a  feasible  solution  when  considering  a  generic  low-cost,  rapidly 
deployable,  and  highly  customizable  optical  signal  processing  framework  for  contributing  SSA  space-based  sensor 
demonstration  projects  under  the  OrbitOutlook  program.  GOSPA  was  designed  around  real-world  sensor  operations 
which  is  why  it  includes  support  for  noisy  data  (low  target  SNR)  and  large  uncertainties  in  sensor  pointing 
knowledge. 

With  only  slight  modifications  to  the  space-based  platform  propagation  routine,  GOSPA  will  support  ground-based 
optical  sensors.  Its  capabilities  can  also  be  adapted  to  process  and  predict  the  performance  of  an  ad-hoc  sensor 
network  based  on  the  previous  satellite  ephemeris  data  provided  by  the  Space  Surveillance  Network  (e.g.  both  IOD 
and  OD).  GOSPA  could  also  be  used  for  sensor  performance  trade  studies  of  ground  and  space-based  optical  sensor 
architectures. 

Ongoing  improvements  are  made  by  AFRL/RV  to  the  feature  point  detection  algorithms  in  GEODETICA  and  to 
CAR-MHF.  Since  the  evaluation  of  CAR-MHF,  another  version  has  been  released  with  updated  filters  and 
smoothers  which  may  allow  for  fewer  and  noisier  observations  in  order  to  satisfy  conditions  for  computing  an  IOD. 
Work  can  also  be  performed  to  provide  support  for  handling  RADAR  and  passive-RF  collectors  within  GOSPA. 

This  would  involve  adding  a  CAR-MHF  feature  to  allow  for  the  processing  of  range  and  range-rate  information 
when  it  is  available.  The  systematic  errors  observed  in  the  CAR-MHF  Performance  vs  Pointing  Angle  Uncertainty 
section  can  be  further  reduced  by  calibrating  satellite  force  models,  using  a  smaller  sensor  detector  IFoVs,  and 
experimenting  with  GEODETICA  streak  end-point  detection  vs  centroiding  routines  in  the  image  processing  portion 
of  GOSPA. 
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